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We study an interquark QQ potential for the charmonium system, that is determined from the 
the equal-time and Coulomb gauge QQ Bethe-Salpeter (BS) wavefunction through the effective 
Schrodinger equation. This novel approach enables us to evaluate a kinetic heavy quark mass 
nig and a proper interquark potential at finite quark mass rag, which receives all orders of 1/rag 
corrections on the static QQ potential from Wilson loops, simultaneously. Precise information of 
the interquark potential for both charmonium and bottomonium states directly from lattice QCD 
provides us a chance to improve quark potential models, where the spin-independent interquark 
potential is phenomenologically described by the Cornell potential and the spin-dependent parts 
are deduced within the framework of perturbative QCD, from first-principles calculations. In 
this study, calculations are carried out in both quenched and dynamical fermion simulations. We 
first demonstrate that the interquark potential at finite quark mass calculated by the BS amplitude 
method smoothly approaches the conventional static heavy quark potential from Wilson loops in 
the infinitely heavy quark limit within quenched lattice QCD simulations. Secondly, we determine 
both spin-independent and -dependent parts of the interquark potential for the charmonium system 
in 2+1 flavor dynamical lattice QCD using the PACS-CS gauge configurations at the lightest pion 
mass, M % = 156 MeV. 
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1. Introduction 

The linearly rising interquark potential in QCD plays an essential role in the formation of 
hadrons. Indeed, the quarkonium states such as charmonia and bottomonia can be described well 
by quark potential models [1, 2, 3], where the Coulomb plus linear potential, the so-called Cornell 
potential, is phenomenologically adopted as the spin-independent central potential. One of the 
major successes of lattice QCD is to demonstrate that the static Wilson loop gives a confining 
potential between infinitely heavy quark (Q) and antiquark (Q), which support the phenomenology 
of confining quark interactions in the heavy QQ system [4]. 

As for the spin-dependent potential, spin-spin, tensor and spin-orbit terms of the interquark 
potential can be identified as relativistic corrections to the static QQ potential, which are classi- 
fied in powers of 1 /nig within a framework called potential non-relativistic QCD (pNRQCD) [5]. 
However, although the spin-spin potential have been precisely calculated as one of the next-leading 
order corrections in quenched lattice QCD [6], their attractive spin-spin potential does not even 
qualitatively agree with the corresponding one in quark potential models, where the repulsive spin- 
spin interaction is phenomenologically required by heavy quarkonium spectroscopy. The spin-spin 
interaction calculated within the Wilson loop formalism seems to yield wrong mass ordering among 
hyperfine multiplets [6]. In potential models, functional forms of the spin-dependent terms are ba- 
sically determined by perturbative one-gluon exchange [2]. Such the Fermi-Breit type potential, 
which appears at the order of 1 / tHq within pNRQCD, would have validity only at short distances 
and also in the vicinity of thq = °°. The l/mg expansion is not formally applicable at the charm 
quark mass. Therefore, properties of higher-mass charmonium states predicated in potential models 
may suffer from large uncertainties in this sense. 

Under these circumstances, we have succeeded to determine a proper interquark potential at 
finite quark mass from the equal-time and Coulomb gauge Bethe-Salpeter (BS) amplitude through 
an effective Schrodinger equation [7, 8]. (See also Ref. [9]). In this proceedings, we will discuss the 
quark mass dependence of the spin-independent interquark potential in quenched lattice QCD [7] 
and then show results of both spin-independent and -dependent parts of the charmonium potential 
in 2+1 flavor full QCD with almost physical quark masses [8]. 

2. Formulation 

Let us briefly review the new method utilized here to calculate the interquark potential with 
the finite quark mass. As the first step, we consider the following equal-time QQ BS wavefunction 
in the Coulomb gauge for the quarkonium states. 

Mr) = £(o|G(x)re(x+r)iee;/ PC ), (2.1) 

X 

where r is the relative coordinate of two quarks and Y is any of the 16 Dirac y matrices. Practically, 
the BS wavefunction can be extracted from the following four-point correlation function 

£ (0\Q(x,t)TQ(x + r,t) {Q(x' ,t s )FQ(y' ,1^ |0> 

x,x',y' 

= ^A n (0\Q(x)rQ(x + r)\n)e- M n^ ^ A Q ^{r)e- M ^\ (2.2) 

x n 
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at the large Euclidean time from source location (t s ). Here both quark and anti-quark fields at t s 
are separately averaged in space as wall sources. M„ denotes a mass of the «-th quarkonium state 
\n) in a given J PC channel. For instance, Y is chosen to be 75 and y, to obtain the rest mass of the 
pseudo-scalar (PS) state (J PC = 0~ + ) and vector (V) state (J PC = ), respectively. 

The BS wavefunction satisfies the Schrodingier equation with a non-local potential U [10] 



V 2 

m Q 



0r(r) + f dt f U(rJ)fr(T f )=E T (k(T), ( 2 - 3 ) 



where ihq denotes the quark kinetic mass. The energy eigenvalue Ey of the stationary Schrodinger 
equation is supposed to be My — 2niQ. If the relative quark velocity v = | V/mg| is small as v <C 1, 
the non-local potential U can generally expand in terms of the velocity v as U(r',r) = {V(r) + 
V s {r)S Q ■ Sq + V T (r)S l2 + V LS (r)L • S + ^(v 2 )}S(r' - r), where S n = (S G • ?)(Sg • r) - S Q ■ Sq/3 
with r = r/r, S = Sq + Sq and L = r x (—/V) [10]. Here, V, Vs, Vt and Vys represent the spin- 
independent central, spin-spin, tensor and spin-orbit potentials, respectively. 

In this study, we focus only on the 5-wave states. Thus we perform an appropriate projection 
with respect to the discrete rotation, which provides the BS wavefunction projected in the 
representation. The Schrodinger equation for the projected BS wavefunction 0r(r) is reduced to 

( - — + V(r) + S G • SqV s (r) W(r) = £ r «/> r (r) (2.4) 

at the leading order of the v-expansion. The spin operator Sq ■ Sq can be easily replaced by ex- 
pectation values. As a result, both spin-independent and -dependent interquark potentials can be 
separately evaluated through a linear combination of Eq.(2.4) calculated for PS and V channels as 

V(r)-E + J_pV 2 Mr) + l V 2 <Mrn 25 



1 ( V 2 <Mr) V 2 0ps(r ) 
m Q \ <j>v(r) (j)ps{r) 



Vs(r) = E hyp + J- \ ^> - \ , (2.6) 



where £" ave = M ave — 2mq and Eh yp = My — Mps. The mass M ave denotes the spin-averaged mass as 
|Mps + \My. The derivative V 2 is defined by the discrete Laplacian with nearest-neighbor points. 

Note here that quark kinetic mass niQ is essentially required in the definition of the interquark 
potentials. Under a simple, but reasonable assumption as lim r ^oo Vs(r) = 0, which implies that 
there is no long-range correlation and no irrelevant constant term in the spin-dependent potential, 
we can obtain the quark kinetic mass from the following formula, 

-1 f V 2 <Mr) V 2 <fr>s(r) 

-hyp 



m Q = lim - — <^ } . (2.7) 

r^°°E hyp { 0y(r) 0ps(r) 



3. Numerical results 

We have performed lattice QCD simulations in both quenched and full QCD in this study. 
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Figure 1: (left) A typical result of V 2 0v/0v — V 2 0ps/0ps as a function of spacial distance r. (right) The 
interquark potential calculated from the QQ BS amplitude at finite quark masses covering the range from 
1.0 to 3.6 GeV. Each curve represents the fit result with the Cornell parametrization. 




Figure 2: The quark-mass dependence of A/ a , A, and a as functions of 1 / mg . We perform the extrapolation 
towards the mg — > °° limit (solid curves) of A/ a, A, and a with a simple polynomial function in 1 / mg. For 
(7, a linear fit with respect to 1 /mg is enough to describe the data with reasonable % 2 /d.o.f, while quadratic 
fits are used for A/ a and A. The results given by Wilson loops are also included as open circles. . 



3.1 Nf = quenched QCD simulation 

In quenched lattice QCD simulations, we use a lattice size of L 3 x T = 32 3 x 48 with the single 
plaquette gauge action at /3 = 6/g 2 = 6.0, which corresponds to a lattice cutoff of aT x f»2.1 GeV. 
The spatial lattice size corresponds to La 3 fm. We fix the lattice to Coulomb gauge. The heavy- 
quark propagators are computed using the relativistic heavy-quark (RHQ) action with relevant 
one-loop coefficients of the RHQ [1 1, 12], which can remove large discretization errors introduced 
by large quark mass. To examine the infinitely heavy-quark limit, we adopt the six values of the 
hopping parameter ff, which cover the range of the spin-averaged mass of 15 quarkonium states 
Mive = \ C^ps + 3My) = 1 .97-5.86 GeV. We calculate quark propagators with a wall source located 
at t s /a = 4. Our results are analyzed on 150 configurations for every hopping parameters. 



4 



Interquark potential for the charmonium system with almost physical quark masses 



Taichi Kawanai 



Table 1: Summary of RHQ parameters calibrated for charmonium system in dynamical QCD simulation. 

*c V Cb c E 

RHQ parameters 0.10819 1.2153 1.2131 2.0268 1.7911 



First, in Fig. 1, we plot a difference of ratios of V 2 0v/0v and V 2 0p S /0p S as a function of 
spatial distance r at K = 0.10190, which is close to the charm quark mass, as a typical example. 
The ratios of V 2 <pr/<)>r are evaluated by a weighted average of data points in the range of (t — 
tsrc)/a = 21 - 23. At a glance, the value of V 2 0v/0v — V 2 0p S /0p S certainly reaches a nonzero 
constant value at large distances, which turns out to be the value of — mgAEhyp- We then obtain the 
quark kinetic masses from the long-distance asymptotic values of V 2 0v/0v — V 2 0p S /0p S divided 
by the measured hyperfine splitting A.Ehyp- 

Using the quark kinetic mass determined here, we can properly calculate the spin-independent 
interquark potential from the BS wavefunctions. Figure 1 displays results of our potential obtained 
at several quark masses. For clarity of the figure, the constant energy shift E ave is not subtracted. 
The resulting QQ potentials at finite quark masses exhibit the linearly rising potential at large 
distances and the Coulomb-like potential at short distances as originally reported in Ref [9]. 

We simply adopt the Cornell parametrization for fitting our data V (r) = —A /r + or + Vo with 
the Coulombic coefficient A, the string tension a, and a constant Vq. All fits are performed over the 
range 1 < r/a < 1 1. In Fig. 2, we show the quark-mass dependence of A/ a, A and a as functions 
of 1/rriQ. First, regardless of the definition of mq, the ratio of A/ a in the top figure indicates that 
the QQ potential calculated from the BS wavefunction smoothly approaches the potential obtained 
from Wilson loops in the infinitely heavy-quark limit. If we pay attention to the quark-mass de- 
pendence of each of the Cornell parameters separately, we observe that, although the Coulombic 
parameter A depends on the quark mass significantly, there is no appreciable dependence of the 
quark mass on the string tension a. Their extrapolated values at tkq — > °° are again consistent with 
those of the Wilson loop result. The extrapolation curve are also displayed as a solid curve in Fig. 2. 

3.2 N f = 2 + 1 full QCD simulation 

Full QCD simulations are also carried out by using 2+ 1 flavor gauge configurations generated 
by PACS-CS Collaboration on lattice of size 32 3 x 64 with the Iwasaki gauge action at j8 = 1.9, 
which corresponds to a comparable lattice cutoff of a~' « 2.2 GeV [13]. The spatial lattice size 
corresponds to L pa 3 fm. The hopping parameters for the light sea quarks {K ut j,K s }= {0.13781, 
0.13640} correspond to M K = 156(7) MeV and M K = 554(2) MeV [13]. Our results are analyzed 
on all 198 gauge configurations. We also employ RHQ action to compute the heavy quark propa- 
gator, which has five parameters K c , v, r s , cb and ce- The parameters r s , cb and ce are determined 
by tadpole improved one-loop perturbation theory [12]. For v, we use a non-perturbatively deter- 
mined value, which is adjusted as c 2 ff = 1 in the dispersion relation E 2 (p 2 ) = M 2 +c 2 ff |p| 2 for the 
spin-averaged 15 charmonium state. We choose K c to reproduce the experimental spin-averaged 
mass of IS charmonium states M^l = 3.0678(3) GeV. As a result, the relevant speed of light, 
c 2 ff = 1.04(5), and the spin-averaged IS charmonium mass, M ave = 3.0702(9) GeV, is observed 
with our RHQ parameters summarized in Table 1. To increase statistics for BS wavefunction, 
we have computed charm quark propagators with two wall sources located at different time slices 
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Figure 3: Spin-independent (left) and spin-dependent (right) parts of the central charmonium potential. The 
solid curves show the phenomenological potentials adopted in a NRp model [3]. The dashed curve appeared 
in the left panel represent a fit result with the Cornell parameterization. 



Table 2: Summary of the Cornell parameters and the quark mass determined from lattice QCD. For com- 
parison, the corresponding values adopted in a non-relativistic potential (NRp) model [3] are also included. 





This work 


Polyakov lines 


NRp model 


A 


0.861(17) 


0.403(24) 


0.7281 


Vo [GeV] 


0.394(7) 


0.462(4) 


0.3775 


m Q [GeV] 


1.74(3) 


DO 


1.4794 



t s /a = 6 and 57 and fold them together to create a single four-point correlation function. Then 
the measured hyperfine mass splitting Mh yp = 0.1137(8) GeV is close to the experimental value 
= 0.1166(12) GeV. 

First, we show a result of the spin-independent charmonium potential V(r) with the fitted curve 
(dashed curve) in Fig. 3, where the constant term is subtracted to set V(ro) = with the Sommer 
scale, ro ~ 0.5 fm. The Cornell parametrization is also simply adopted for fit. We have earned 
out correlated % 2 fits with full covariance matrix for on-axis data over range 4 < r/a < 10. The 
fitting results are listed in Table 2. The quoted errors represent only the statistical errors given by 
the jack-knife analysis. The phenomenological potential used in NRp models [3] is also plotted as 
a solid curve for comparison in Fig. 3. Although the charmonium potential obtained from lattice 
QCD is quite similar to the one in the NRp models, the string tension of the charmonium potential 
is slightly stronger than the phenomenological one. Therefore our result indicates that there are 
only minor modifications required for the spin-independent central potential in the NRp models. 
Moreover, it seems that a large gap for the Coulombic coefficients between the conventional static 
potential and the phenomenological potential is filled by our new approach, where all orders of 
1 /rriQ corrections are nonperturbatively accounted for. 

In Fig. 3, we next show the spin-spin term of the charmonium potential and the corresponding 
phenomenological one found in Ref. [3]. Our spin-spin potential exhibits the short range repulsive 
interaction, which is required by the charmonium spectroscopy. It should be reminded that the 
Wilson loop approach fails to reproduce the correct behavior of the spin-spin interaction, since the 
leading-order spin-spin potential classified in pNRQCD becomes attractive at short distances [6]. 
In contrast of the case of the spin-independent potential, the spin-spin potential obtained here is 
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absolutely different from a repulsive 8-f unction potential generated by perturbative one-gluon ex- 
change. Indeed, the finite-range spin-spin potential described by the Gaussian form is adopted in 
Ref. [3], where many properties of conventional charmonium states at higher masses are predicted. 
This phenomenological spin-spin potential is also plotted in Fig. 3 for comparison. There still 
remains a slight difference between the spin-spin potential from first principles QCD and the phe- 
nomenological one. In this sense, the reliable spin-dependent potential derived from lattice QCD 
can provide new and valuable information to the NRp models. 

4. Summary 

We have proposed the new method to determine the interquark potential at finite quark mass 
from lattice QCD. Using quenched lattice QCD, we have demonstrated that the spin-independent 
central potential defined in this method smoothly approaches the static QQ potential given by Wil- 
son loops in the infinitely heavy-quark limit. In dynamical lattice QCD simulations, we have stud- 
ied both spin-independent and -dependent parts of the charmonium potential. The spin-independent 
charmonium potential obtained from lattice QCD with almost physical quark masses is quite sim- 
ilar to the one used in the NRp models. The spin-spin potential properly exhibits the short range 
repulsive interaction. Its r-dependence, however, is slightly different from the phenomenological 
one adopted in Ref. [3], Therefore, our charmonium potential derived from first principles QCD 
suggest that properties of higher-mass charmonium states predicted in the NRp models may change. 

We acknowledge the PACS-CS collaboration and ILDG/JLDG for providing us with the gauge 
configurations. We would also like to thank H. Iida, Y. Ikeda and T. Hatsuda for fruitful discussions. 
This work was partially supported by JSPS/MEXT Grants-in-Aid (No. 22-7653, No. 19540265, 
No. 21105504 and No. 23540284). 
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